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COMPUTER EXPERIMENTS ON THE MICROSCOPIC BEHAVIOR 


OF ONE -DIMENSIONAL PLASMAS* 

By Richard H. Weinstein and Belinda H. Adams 
Langley Research Center 

SUMMARY 

The detailed motions of a system of 1000 particles simulating a one-dimensional 
electrostatic plasma have been followed by using the IBM 7040/7094 Direct Coupled 
System at the Langley Research Center. The particle interactions are collective and, 
as the calculation proceeds stepwise through time, the position-velocity histories of all 
the particles completely specify the successive states of the system. These data consti- 
tute an "experiment, " upon which measurements can be made. The velocity distribution 
function, the Fourier harmonics of the potential energy, the spatial correlations of the 
electric field, and the drag and velocity diffusion of "test particles" are measured for the 
Maxwellian and uniform velocity distributions. Comparison of the present work with both 
earlier experiments and the predictions of zeroeth and first-order plasma kinetic theory 
have yielded generally good agreement for all measurements in both equilibrium and non- 
equilibrium plasmas. The detailed behavior of the long wavelength harmonics of the 
potential energy is, however, unresolved for regions where the assumptions of the 
zeroeth-order theory are not satisfied. 

INTRODUCTION 

The large information- storage capacities and high calculating speeds of the digital 
computer have made possible a type of "computer experiment" which can provide rigidly 
controllable simulation of N-body interactions such as occur in plasmas. In contrast to 
Monte- Carlo techniques which are at best iterative, the computer experiments are 
dynamically self-consistent because they treat the simultaneous motion of all N parti- 
cles due to the instantaneous, self-consistent Coulomb field of the charged particles. 

This method can be viewed as an experiment against which the predictions of theory may 
be checked. 


*The information presented herein was offered as a thesis in partial fulfillment of 
the requirements for the degree of Master of Science in Physics, Virginia Polytechnic 
Institute, Blacksburg, Virginia, June 1967. 


Although some three-dimensional problems (ref. 1) have been studied, the number 
of particles treated is generally modest, both because of the large amount of information 
storage required for each particle and because of the complexity of the calculations in 
several dimensions (for example, solution of the self-consistent potential). The sim- 
plicity of one-dimensional models allows the use of large numbers of particles and 
thereby improves the statistical accuracy of the results. A one-dimensional model can 
be constructed by visualizing the particles as parallel sheets and allowing them to pass 
through (that is, cross) each other. Such models, used in the one- dimensional plasma 
diode experiments initiated by Hartree and Nicholson (1943, unpublished) and continued 
by Lomax (ref. 2), Birdsall and Bridges (refs. 3 and 4), and by Burger (ref. 5), can pro- 
vide realistic dynamic simulation of a laboratory experiment and thereby serve as a 
guide to developing a theory. Conversely, the theoretically predicted two- stream insta- 
bility was first observed ’’experimentally’' by Buneman (ref. 6) in a computer experiment. 
Monte-Carlo calculations, such as those of Goldstein (ref. 7) solve only for the steady 
state of the diode, and self-consistency is approached iteratively. Quasi-one-dimensional 
models have been used for problems involving magnetic fields, such as the ion- cyclotron 
wave calculation of Hasegawa and Birdsall (ref. 8). A rather sophisticated one- 
dimensional model treated by Dawson (ref. 9) has included self-consistent electromagnetic 
radiation. Several two-dimensional experiments have been used as engineering tools; for 
example, Buneman and Kooyers (ref. 10) studied the ion engine neutralization problem, and 
Yu, Kooyers, and Buneman (ref. 11) used a model of the magnetron which could, in princi- 
ple, be used to design tubes. A complete summary of recent work in both computer experi- 
ments and simulation of plasmas by computer solution of the governing equations can be 
found in reference 12. 

This work follows the approaches of Eldridge and Peix (refs. 13 and 14) and Dawson 
(refs. 15 and 16), who were concerned with checking predictions of the Vlasov-Landau 
theory, rather than simulating an experimental situation. An approximate computer pro- 
gram is used in which several particle crossings are allowed per time step of the calcu- 
lations; this approach has made it practical to treat larger numbers of particles per 
Debye distance (nD), that is, more nearly collisionless (Vlasov) plasmas, and to follow 
the plasmas for long times. The accuracy of the calculation is checked both through 
conservation of energy and by its agreement with the Vlasov theory predictions checked 
in earlier experiments with exact models (refs. 13 to 16). 

Studies of nonequilibrium plasmas are the principal extension of this work over that 
of references 13 to 16. Eldridge and Feix (ref. 14) predicted a stationary state or 
"metaequilibrium" in one dimension for any stable distribution function; this phenomenon 
was first observed by Dawson (ref. 16). The larger values of nD used in this report 
show the expected improvement in stationarity of the distribution function. 
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A brief review of the kinetic theory is provided as an introduction to the potential 
energy harmonic and spatial correlation theory. Measurements of the Fourier harmonics 
of the potential energy and the electric-field spatial correlation are repeated for a non- 
equilibrium plasma. The first-order kinetic theory is presented next, although the meta- 
equilibrium, a first-order result, underlies even the Vlasov nonequilibrium studies. The 
drag and velocity diffusion of test particles are measured both to check the scaling pre- 
dictions of the kinetic theory and its quantitative accuracy. 

SYMBOLS 

A normalizing coefficient of one-particle distribution function 

a limiting velocity of square distribution 

C autocorrelation function 

c mode number of Fourier component 

D Debye distance 

D 2^3 Debye distance analog for f2(v) and f3(v), respectively 

electric field 

F 

dimensionless electric field, - — 

47TO- 

N-particle distribution function 
F reduced many-particle distribution function 

f one -particle distribution function 

g graininess parameter, ^ 

Im( ) imaginary part of variable 

k wave number 

L length of system 
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contours shown in figure 4 


^P’^oo 

m mass per unit area of sheet particles 

m e mass per unit area of electron sheets 

N total number of particles 

n system particle density, N/L 

P pair distribution function 

q number of particles described by distribution function Fq 

Re( ) real part of variable 


r 

s 

T 

t 

V 

v 

V T 

W 

x 

y 

Oi 


exponential power index of one-particle distribution function 
Laplace transform variable, complex frequency 
dimensionless time, Wpt 
time 



thermal velocity or root-mean- square velocity 


energy 

position 

dummy variable 
acceleration 
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/3 normalization factor in one-particle distribution function given in dimensions 



y order of expansion parameter 

e plasma dielectric constant 

?7 total energy growth rate 

6 angle in complex plane 

X electric-field correlation length 

I u charge-to-mass ratio of particles normalized to electron values, 

p charge number density 

ct absolute value of the charge per unit area of electron sheets 

Wvjt 

r nondimensional relaxation time, — £- 

nD 

<p Fourier- Laplace transform of one-particle distribution function 

X Fokker- Planck drag coefficient 

Fokker- Planck diffusion coefficient 

radian frequency 

4„ ncr 2 

plasma frequency, — ^ — 

Subscripts: 

h,i,j indexing parameters 

A bar over a symbol denotes a nondimensional value; an arrow over a symbol 
denotes a vector. Average values are denoted by < >. Primes denote perturbed quan- 
tities and an asterisk denotes a complex conjugate. 


Q, = 

Cl) 


= ^p 
CD k 


Cl)v 
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EXPERIMENTAL SYSTEM 


Model and Equations of Motion 

The plasma is simulated by following the motions of large, plane, parallel sheets 
confined to one-dimensional motions along an axis normal to the sheets. Equal numbers 
of positively and negatively charged sheets are used so that the total system is electri- 
cally neutral. The total number of sheets in the system is taken to be N and all sheets 
have the same mass density. Each sheet is acted on by the collective coulomb forces 
of all the other sheets. The electric field in the system is a step function; it is constant 
between sheets and changes by ±4 ttct in crossing a sheet, where cr is the charge per 
unit area of an electron sheet. A field configuration for a typical plasma having eight 
particles is shown in figure 1(a), where the circles represent the value of the field at the 
position of the sheets. 

The sheets, or "particles" as they will be called henceforth, are allowed to pass 
through each other (that is, cross). The system is considered to be periodic; thus, the 
motion of the particles is followed only between two imaginary boundaries at x = 0 and 
x = L. This condition is physically equivalent to representing an infinite plasma by a 
periodic structure of fundamental wavelength L, and thereby placing a maximum wave- 
length constraint on the Fourier representation of system quantities, such as density and 
potential energy. To describe an infinite plasma with this model successfully, it is 
required that L/D » 1, where D is the Debye distance, the characteristic scale length 
for a plasma. For computational purposes, the system may also be viewed as having its 
two ends (x = 0 and x = L) tied together. 

Particles leaving this interval are accounted for by cumulative charge sums kept 
on the wall crossed where erg is the charge at x = 0 and cr^ +1 is the charge at 
x = L. When a particle crosses a boundary, its "twin" (a particle with the same charge 
and velocity vector) is injected at the opposite wall along with a compensating charge 
kept on the opposite wall. This condition is shown in figure 1(b), where particle 8 has 
crossed the boundary at x = L and has been reinserted at x = A. The field change 
across erg is - 47 rcr; thus, +4tto and -47rcr are added, respectively, to erg and cr N+ i, 
and the twin particle 8' appears near x = 0. This manipulation of boundary charges 
maintains the self-consistency of the field in the plasma regardless of events at the 
boundary. Thus, for example, the field seen by particle 7 remains at +2ir<y; without the 
addition of Og and Ojq. j, this field would have been -2iro after the crossing, and the 
conservation of electrostatic (that is, potential) energy in the system would have been 
destroyed. 

Because of charge neutrality, the electric field has the same value at x = 0 and 
x = L, although it may be nonzero because of fluctuations of erg and cr N+1 . At a given 
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time t, the force per unit area experienced by a representative sheet at position xj is 

m *j = a j E (V) (1) 

where ctj is the charge per unit area, E^Xj,t^ is the electric field, m is the mass 
per unit area, and xj is the acceleration of the jth sheet. The self-consistent field 
seen by the sheet is given by a sum over all the other sheets; that is 


E 


(V) = 


2tt 


N+l 

^ 0^ sgn(xi(t) - xj(t)) 
i=0 


(2a) 


where 


sgn(x) = 


(x>0) 
(x = 0) 
(x < 0) 


The charge densities are henceforth referred to as charges. Because of charge 
neutrality in the interval 0 to L, E^Xj,t^ may be more simply calculated as 


E^Xj,tj = 477 cr^ + 27r<7j (2b) 

Xi<Xj 

where the summation extends over all i such that x^ < Xj. 

The field seen by a particle changes, in general, when it crosses a neighbor, as is 
shown in figure 1(c). Here, particles 2 and 3 have crossed, and the field seen by parti- 
cle 2 has changed from -2770' to -6770. To calculate the exact evolution of the system, 
it is necessary to calculate the time at which this crossing occurs, to advance the parti- 
cles concerned to the crossing point, and then to recalculate the electric field to be used 
for the next step of the motion. Because the number of crossings that must be considered 
increases both with the total number of particles and nD, the range of problems feasible 
with such an exact program is somewhat limited because of the long calculation times 
necessary. Such an exact program was used by Eldridge and Feix (refs. 13 and 14). 


This paper uses an approximate program based on the integral of equation (1) over 
a small fixed time step At; the initial field E(Xj is taken as the forcing term for the 
motion over the time step and the field is recalculated only at the end of the time step. 

If the size of the time step is small compared with the inverse plasma frequency, which 
serves as a fundamental scaling time for the plasma, few crossings leading to only small 
errors in the particle motion will have occurred. These crossings destroy the strict 
conservation of energy possible with the exact program within computer truncation error, 
and the extent of departure from strict conservation provides a criterion of acceptability 
for the approximate program. The errors in particle motion are due to the fact that 
crossings make the electric field a function of time during the motion. The electric field 
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E^xj,tj will, in fact be discontinuous, but if the particle motion is to be recalculated only 
at the end of each time step, one may smooth the field over a time At and obtain cor- 
rection terms from a Taylor series expansion; for example, with At = t 2 - t j, 

E (x.,t 2 j = + E(xj,tj)At + | E^tjjAt 2 + . . . 


These time derivatives may be approximated by their backward difference equivalents, 
for example, 



E ( x i ,t l) “ E ( x j jt l " At ) 

At 


(3) 


Integration of equation (1) from t = 0 to t = At gives for the position and velocity of 
the jth particle, 


'i( l l + 4t ) = x j(‘l) + Vj(ti)it + § 5 E(xj,‘l)it 2 + i E(xj,t x )it 3 + . . . (4a) 


v i(‘l + 4t ) = T j(‘l) + | %‘lK + . . j (4b) 


The time increment At is taken to be small and terms of order At 2 and higher 
are neglected; as a result, the equations of motion can be truncated after the (At) 2 term. 
To this order, the time dependence of E(xj,t) gives rise to the (At) 2 term in the veloc- 
ity equation, which would not appear in a constant- force situation; the strong effect of this 
term on the conservation of energy is shown. 


Equations (4a) and (4b) can be simplified by normalization in terms of characteris- 
tic parameters of the plasma, namely, D the Debye distance, Vq> the thermal velocity 

(the root mean square velocity for an equilibrium plasma), and o>p s the plasma 

Vt> 

frequency. These quantities are interrelated by the definition D = ■—% Also of interest 

P 

are the following dimensionless quantities: 



With the normalizations, 


AT = Atwp 



E ( x j> t ) 

47TCT 
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and 



equations (4a) and (4b) become 


D AT + | -ji E(xj,T AT 2 

(5a) 

j’T^AT + — E^Xj,T^AT^j 

(5b) 


The position x^ is not normalized to the Debye distance D for reasons of convenience 

J "NT 

and' in the computer program the density n = — is taken to be unity. 


Initial Condition 

The initial state of the plasma is one of uniform particle density; the particles are 
equidistant from one another in the interval x = 0 to x = L and electrons and ions 
alternate. This configuration is the minimum potential energy state of the uniform system 
and was chosen because of its computational convenience. Individual particle velocities 
are initially chosen from a pseudo- random number generator. Both velocity distribu- 
tions used in the experiments described here are of the form f(v) = A exp(-/3v r ), where 
r is a positive number. The quantities A and /3 are independent of v but are deter- 
mined by the normalization of f(v) and the second moment of f(v), respectively. The 
values r = 2 and r = °° give the Maxwellian [fi(v)j and uniform ^3 (v)j distributions, 
respectively. These velocity distributions are shown in figure 2, along with f2(v), the 
Druyvesteyn distribution (r = 4) which proves to be an interesting intermediate step in the 
transition of the uniform distribution f3(v) to equilibrium. 

All distributions are normalized to have the same total kinetic energy for a given 
value of D. This requirement dictates the sharp cutoff of the uniform distribution 
yielding f3(v) which will henceforth be referred to as the "square" distribution. The 
Maxwellian distribution has an effective cutoff less than 2.5vq' because of the finite num- 
ber of particles considered and the exponentially decreasing form of the distribution 
function. 

It should be noted that for each type of velocity distribution (that is, Maxwellian or 
square), the same set of pseudo- random numbers is used for all values of nD. Changing 
nD is, however, sufficient to give different time histories to each plasma (At being held 
constant) because, for example, a particle traveling with the thermal velocity through a 
plasma of density n crosses nD particles per AT. Thus the electric field seen by 
each particle is, in general, different even after the first time step, because the different 
numbers of crossings for each value of nD give different charge configurations. The 
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scaling of first-order plasma kinetics with nD is not, therefore, a trivial result of 
scaling the equations of motion, but is shown to be a prediction of plasma kinetic theory. 

Computational Information 

The computer time necessary to advance all particles through one cycle (that is, to 
update the system from time T to time T + AT) depends on both the system size and 
the plasma parameters. For the range of conditions studied here, this time is empiri- 
cally given for the IBM 7040/7094 Directly Coupled System in seconds as 

T cycle = 5N x 10- 4 + 0.89NnD AT x 10‘4 

The first term in this equation represents the time necessary to calculate the new posi- 
tions and velocities for all N particles from stored or calculated quantities. The sec- 
ond term is due to calculation of the electric field. Because the field seen by a particle 
depends on the excess charge on either side of it, a table of identification numbers giving 
the order of particles is recomputed each time the particle positions are recalculated. 

The time necessary to accomplish this computation obviously depends on the size 
of the time step (longer times yield more crossings), and the necessity of scaling with 
nD can be seen from the fact that a particle traveling at the thermal velocity crosses 
(on the average) nD particles in a time Typically, for a plasma of 1000 particles 

with nD = 40 and AT = 0.1, the cycle time is =0.83 second. This time represents con- 
siderable improvement over an exact program, which has a speed comparable to the 
approximate program for nD = 10, but the time unfortunately grows as ND rather than 
as nD. As the calculations proceed stepwise through time, the position and velocity of 
each particle along with the total kinetic and potential energy are stored on magnetic tape 
at each cycle or desired multiple thereof. These data are used as inputs for the mea- 
surements described subsequently. Some of the techniques used in the plasma generation 
and measurement programs are discussed in the appendix. 

Conservation of Energy 

The energy calculation shows the effect of including E in equation (5b). Whereas 
the conservation of energy for the exact program is precise (within the limits of com- 
puter precision), the total energy grows slowly as a function of time with the approxi- 
mate program.. The rate of total energy increase, scaled to the initial energy and the 
plasma frequency, is taken as 

. ( AW tof-p 

" CMo 

where ^W^-q^q is the total initial energy, AWtot is the difference between final and 
initial energy. 
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Figure 3 shows this growth rate as a function of the size of the time step for a 
range of values of nD. Note that different values of nD correspond to plasmas at dif- 
ferent temperatures, for = f^\ cc (Kinetic energy). If the motion in equation (5b) is 

V^P / 

calculated without E, the empirically determined form of the energy growth rate is found 
to be 


V = 


AT 

2(nD - 1) 


This expression is valid over the range of AT shown in figure 3. Comparison between 
curves with and without E for corresponding values of nD shows that one can increase 
the size of AT by about a factor of four for a specified accuracy; this increase in AT 
proportionately decreases the time needed to advance the system one time step. This 
decrease is the justification for using an approximate program for Vlasov plasmas. 
Although the number of crossings and, therefore, the deviations from strict accuracy 
increase with nD, the smoothing of fluctuations in the fluid limit aids the gross physical 
consistency. 


ZEROETH-ORDER THEORY AND MEASUREMENTS 
The Plasma Kinetic Theory 

Experiments using the plasma model described have been compared with predic- 
tions of both zeroeth and first-order kinetic theory derived from the Bogoliubov-Born- 
Green-Kirkwood-Yvonne (BBGKY) heirarchy (ref. 17) outlined briefly below for a one- 
dimensional system (ref. 18). The starting point is the Liouville theorem, which 
describes the motion of an ensemble of identical systems in 2N- dimensional phase space 
(that is, x,v space), where N is the number of particles in each system. The state 
of each system in the ensemble is denoted by one point in phase space. If 

<r -' N 

( x l> v l> • • •» X N> v N,t). n , ^i dv i 

represents the probability of finding a system at a particular point in phase space at a 
time t, with particle 1 between Xj and x^ + dxj, and so forth, and is normalized so 
that 

p°° poo Cp N 

\ . . . \ ~Jt II dx4 dv^ = 1 

CO O_oo i = l A 1 

then the Liouville theorem for the system states that 
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( 6 ) 



where is the acceleration experienced by the ith particle. 

If one defines the reduced q-particle distribution function as 

pOO pOO fy, N 

F a = L q \ . . . \ n dXi dVi (7) 

H J-co J-oo i=q+i 

where L is the length of the system, then multiplication of equation (6) by 

_ N 

L q n dxj dVi 
i=q+l 


and integration over N - q coordinates yields a linked set of equations giving Fq in 
terms of F^ + j; that is, use of the indistinguishability of the N - q particles and neglect 
of surface integrals at the phase- space boundary yields 




8F 


q+1 


&v< 


( 8 ) 


In equation (8), 2irOj sgn^Xj - Xjj is the acceleration of the ith parti- 

cle due to its interaction with the jth particle; that is, external fields are taken to be zero. 
Dependence of Fq on F^ + j arises from a sum over the two-particle interaction poten- 
tial a. „ . i so that in principle, all the higher order functions must be known in order to 


i,q+l 


determine the kinetic equation for any given function Fq. The system must be large 
enough (N L — °°) to insure that the graininess expansion parameter (to be intro- 

duced next) is still large compared with perturbations n’/n where n = ^ = Constant is 
the equilibrium density. 

The graininess expansion is introduced by proceeding to a fluid limit description of 
the plasma. For a large system, one conceptually splits each particle so that the elec- 
tronic charge and mass densities (o and m) become very small, and simultaneously 
defines the limits that cr/m is constant and that nm and ncr remain finite. Although 
this fluid limit is a purely mathematical construction, a real plasma exhibits the same 
type of fluid properties because long-range coulomb forces result in collective behavior. 
Equation (8) may be made dimensionless by introducing the plasma frequency, thermal 
velocity, electron charge, and charge- mass ratio 
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where m e and cr refer to the absolute values of electron quantities. Electric fields 
are normalized in terms of the sheet charge density, that is, E is normalized as 
E = This normalization gives rise to the so-called graininess parameter 


« - (-^f 2 

mmvrj.'y 


_ 1 
nvrp nD 


which, for large nD, may be used as an expansion parameter for the linked equations. 
If the q-particle distribution function Fq is expanded in orders of g, that is, 


F q = 


+ gF q + g 2 F q ^ + . . . g r F q y + 


(9) 


Then substitution into equation (8) yields, after separation by powers of g, dimensionless 
kinetic equations to any desired order y; that is, 


0F^ 

— a^ + 

at 


V 8F q^ 1 r°° 9F 1+1 

Z [i -*T + 2 Vl **+1 % + l Sgn ( X i ' Vl)-^T 


IV V— aF a r ^ 

= - 2 Z ^ Z ff J sgn ( Xi ' x i) 

i=l 3=1 


( 10 ) 


Complications in equation (10) arise because F^ is given both in terms of F^ y_ ^ 

and F q^.j* The ordering in powers of g is equivalent to introducing a hierarchy of 
time scales in which the Xth-order kinetic equation describes 9F q j^ T y where 

r = is the yth-order time scale. As a final step, some relation between F n and 

y C 0 p 4 

Fq + j must be specified. The Mayer cluster expansion (ref. 17) is the usual choice and 

gives the following form: 

F l^ = F 1°^ 1 ) 

F 2 (l,2) = F^(1)F^(2) + P(l,2) 

F 3 (l,2,3) = F^(1)F^(2 )fJ 0) (3) + F^(1)P(2,3) + F^(2)P(1,3) + fJ 0) (3)P(1,2) + T(l,2,3) 

In this expansion, P(l,2) = gF^'(l,2) is the pair correlation function and 

T(l,2,3) = g2Fg2(i > 2,3) is the triplet correlation function. The chain of equations may 
then be formally broken by neglecting small quantities — in this case, any term of order 
g2 or higher. 
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To zeroeth order, one obtains the dimensionless Vlasov equation (ref. 18), 

]f(l) = 0 


d d —(0) a 

— +V., + Mi E -T— 

3t 1 3xi 1 3vi 


(id 


where f (1) = F^(l) is the usual one-particle distribution function and E'®' is the 
self-consistent field determined by Poisson's equation; that is, 

— (0) if 00 

E' = i J dXg dv2^ s ^ n ( x l “ x 2jl(2)- If a homogeneous system is considered, the 

first-order equation for the total one-particle distribution function is, from equation (10) 
with q = 1, y = 1, 


® 02 13X2 dV2 ‘ X2 ) p(1 ’ 2) 


( 12 ) 


and P(l,2), in turn, satisfies the first-order equation for q = 2, y = 1: 


8 8 8 

— + V 1 + V 


at 1 8Xi 2 9X, 


r ? 3 *3 dv 3 s * I f2 ‘ X 3 ) P( 2 ’ 3 >^r| = f - x 2)' 1 l a 2 f (2)^ 


^ p (l,2) + !£ a 3 dx 3 J ^ dv 3 f S gn(x 1 - x 3 )P(l,3)Mll 




(13) 


Subscripts and numbers in parentheses in equation (12) correspond to specific particles. 
The Vlasov equation (eq. (11)) provides the basis for interpreting the potential-energy 
spectrum and electric-field correlation measurements; the first-order kinetic equations 
(eqs. (12) and (13)) provide the relaxation theory. 


Vlasov Theory 


Solution of the Vlasov equation for a homogeneous one-dimensional system (refs. 19 
and 20) is conveniently handled by using the symmetric Fourier- Laplace transform of the 
one-particle distribution function; that is 


<M k,v,s) 


OO °0 

— \ dx \ dt f(x,v,t) exp(ikx) exp(-st) 
2tt -Loo Jq 


and similar transforms for other functions of x and t. For a system composed of 
point particles, the 6-function density in position space is changed by a Fourier transform 
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into a sum of continuous functions in k- space. The Laplace transform on time introduces 
initial conditions into equation (11); that is, if 0(k,v,s) is the Laplace transform of a 
function f(x,v,t), then 

C dx exp(ikx) f dt exp(-st) = s0(k,v,s) -f(k,v,t=0) 

J _ CO Jo ^ 

If one takes 


J oo pc 


dv f(x,v,t) 


1 


one can consider a perturbed one-particle distribution function of the form 

$ = 4>°(v) + $'(k,v,s) 


(14) 


where is the homogeneous velocity distribution function and c£) T (k,v,s) is a small 

perturbation. The integral of this perturbation over phase space must be zero. The 
perturbed distribution function 0 T (k,v,s) gives rise to a perturbed electric field E f (k,s) 
so that substitution of equation (14) into the dimensional Vlasov equation and linearization 


(that is, neglect of the second-order of E T termj gives for ions (+) and electrons (-), 
respectively, 


cj) (k,v,s) = — 

± s - lkv 


0±(k,v,t=O) - (^) E'(k,s) 


9f°(v) 


dV 


(15) 


It is to be noted that f^(v) = <fi®(v) and is the unperturbed distribution function which is a 
function only of velocity. Here the externally applied field is still assumed to be zero and 
the perturbed linearized field is given by the Fourier- Laplace transform of Poisson's 
equation 


J OO ago . v poo poo poo 

dx \ dt exp(ikx - st) ' x? = 4 ttou \ dx \ dt exp(ikx - st) \ (f' - f'jdv 

- oo J — oO J _ GO J _ oO J _ go ' V 


or 

p oo 

-ikE’(k,s) = 47ron J dv|^/ (k,v,s) - <£/_(k,v,s)j = 4irp(k,s) (16) 

where a is the charge density of an electron sheet, and all fields are taken to be zero 
at infinity. 

If both species have the same mass density and distribution function, then combina- 
tion of the electron and ion equations through the use of Poisson's equation gives for 
p(k,s) the net charge density, 
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p(k,s) = 


C£ 

J_oo f 


e(k,s) J -oo s - ikv 


where 


e(k,s) = 1 + 


r°° af%v 

'-00 ^ - iv 
k 


is the plasma dielectric function, which contains a complete description of the collective 
aspects of the plasma (refs. 19 and 20). The Laplace transform variable s is complex, 
and Re(s) represents growth or decay of oscillation amplitudes and Im(s) represents 
a wave propagation frequency co. For Re(s) = 0, e(k,s) = 0 gives a dispersion rela- 
tion for the usual plasma oscillations; recovery of a small Re(s) gives the Landau 
damping (ref. 17). 

For our purposes, it will be more convenient to work with frequency a>, rather 
than with the Laplace variable s. If s = iw, then the Fourier density function can be 
shown to be (refs. 21 and 22) 




^ <°(S) 


:(k,w)| < 


where f^0 is the value of the unperturbed one-particle distribution function at 

v = and co is a real frequency. This relation is very general, as pointed out by 
Rostoker (ref. 21), and represents an extension of the Nyquist fluctuation dissipation 
theorem through the generalized dielectric function e(k,o>). As such, it relates the 
fluctuations to the correlations wherever a stationary distribution function can be defined 
(ref. 22). 

It can be shown (ref. 21) that the Fourier transform of the spectral density function 
is the autocorrelation function C(X,r). If the symmetric Fourier transforms are defined 
as 

E(x,t) = y dk y du> E(k,co) exp£-i(kx + cot)J 

E(x+X,t+r) - dk' J dco’ E(k',u>’) exp£i(k'x + w’t)j exp£-i(k'X) + w'r) j 


use of Poisson's equation (eq. (16)) gives for the autocorrelation function, 
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C (A,t) = <E(x,t) E*(x+A,t+r) > x t = ^ dk J dto <|E(k,co)| 2 > e i ( k ^ +WT ) 

. 8, P dk f “ do, <|p(k ’" >|2> e‘< kx+ " T ) < 20a > 

where the Fourier density functions are given by their ensemble averages. 

For x = t = 0, C(A,r) = T ) XL is the potential energy of the system normal - 


87rmv-p 

ized to twice the kinetic energy; that is, 

C (0>0) „ C(0,0) XL , L <1 E l 2> x,t _ L 


87rmv n 


mv n 


87T 


47tD^ ' j -°o 




OO < 

dco — 

- OO 


P(k,co)|‘ 

no 2 k 2 


Thus, 


C(0,0) 



dk 


<|p(k)| 2 > 

na 2 k^ 



dk W(k) 


(20b) 


In equation (20b) the integration on co has been performed and the remaining integrand 
W(k) gives the Fourier harmonics of the normalized potential energy. If X * 0, one 
obtains C(A,0) the normalized spatial correlation of the electric field, which is related 
to the Debye screening; that is, 


C(X,0) = ± — < E(x,t) E*(x+A,t) > = -±- P dk e ikX (20c) 

87rmVrp 2 47 tD^ -°° na^k* 4 

The quantity E(x) E*(x+A) can be directly measured in the plasma, whereas its 
analytic form can be determined by integration of equations (20b) and (20c) for a specific 
f(v). The distribution function appears both explicitly in the factor and implicitly 


in the e 


: + (k,w) 


factor of < 


|p(k,co)| ' 


>. 


The Fourier Density Function Theory 

For fj(v) and f 2 (v), contour integration in the complex co -plane is the most con- 
venient way of performing the integrations indicated in equations (20b) and (20c), but the 
dielectric function must be more precisely defined. The function e(k,co) can be taken 
to be analytic in either the upper or lower half co -plane but has a discontinuity along the 
entire real axis. For the distribution functions being considered, the only resonances in 
p(k,co) come from zeros of e (k,co), which occur in the upper half co-plane. Therefore, 
the dielectric function was chosen to be analytic in the lower half plane (denoted as 
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e + (k,o>), where the + sign is carried over from the positive real s half-plane) with ana- 
lytic continuation into the upper half u> -plane. The function e + (k,co) is, therefore, an 
entire function of to when defined as follows: 


e + (k,o>) = 1 


e + (k,o>) = 1 


e + (k,co) = 1 


+ ^f dv ff^vj/sv 

k 2 J-oo « - V 


(Im co < 0) 

(21a) 

i w p 2 p r 00 dv 9f0 ( v )/ 8v i7r w p 2 af0 (v) 


(Im co = 0") 

(21b) 

k 2 P J.„ dV <■>.„ +w k 2 dv 

k J 

< 

II 



(Im w > 0) 

(21c) 


k 


In equation (21), Imai = 0" implies the value on the real axis approached from below. 

By using the Plemelj formulas, e + (k,co) is seen to be continuous in the entire co-plane. 

2 

Some simplification in evaluating <|p(k)| > can be accomplished by explicitly sub- 
stituting f®(v) in < | p (k, oj ) | ^>; that is, with f®(v) = Ae"^ r , then 


f0/w\ _ _ 1 8f°(v) 

W rj3v r-1 87 


Jv^ 
k 


If it is noted from equation (21b) that — ^ 


9v 


V= F ffW p 2 


Im e(k,co), then 


w a >-s: 


dco <|p(k,co)| 2 > = -ncr 2 — — 


irr/3o) p 2 


-°° co r 'l |e(k,co)| 2 


= na 2_j^_r 

fr 2 * 


dco 


7rrj3co„ 2 J-oo co 


r-1 


Im 


e(k,o>) 


= na 2 — — 


7rr/3cu 


P 


■s 


*■> -J_ lm(i) (22) 


-oo CO 


where the last step follows because any real convergent integral equals its own principal 
part. This relation may be rewritten as 


< 



> = na 2 


k r 

7rr/3cOp 2 


Im P 



1 

1 

o> r -l 

6_(k,co) 


dco = n 


a 2 


kf 9 Im I(k) 

irrpojp* 


The integral I(k) may then be evaluated by contour integration along the path shown in 
figure 4. The path proceeds along the real co-axis dented along l p into the lower half 
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co-plane to avoid the pole of I(k) at co = 0 and is closed in the lower half-plane along 
l w where |co| — °°. The integrand is analytic within this contour so the Cauchy residue 
theorem gives 



dco 


1 

r i 

C0 r- 1 

e_(k,co)_ 



1 1 C dco 1 

a, 11 - 1 [e(k,co)] co r_ 1 [e(k,co)] 


(23) 


The first case to be treated is the Maxwellian distribution, for which r = 2 and 
there is a simple pole at co = 0. The integral along l «, does not vanish in this case but 
is finite because lim e(k,co) ~ 1, with deviations of order (ref. 21). Thus, with 

jcoj— co Cl/ 

co = p exp id, 


I(k) - -iri 


e.(k,co) 


co=0 0 


p-77 

Jo 1 


i d0 = -7ri 


- 1 


COy 


1 + 


Using the second moment of the distribution function to show that /3i = — for fi(v), 

2v 2 

Vrp 

together with D = — - , then 

COp 




2 . _ _2 k 2 P 2 

^ ^ 9 9 

1 + k z D 2 


■°(v) = f i(v)) 


(24) 


For the Druyvesteyn distribution (r = 4), the integral along loo vanishes because 
of the co -3 dependence of the integrand. Considerable care must, however, be taken 
along the dented contour l p because the contribution from a partially included singular- 
ity is generally finite only for a simple pole. The contribution will be finite in this case 
for several reasons shown subsequently. 


In evaluating I(k,co) along l p , where I(k,co) = - y the analytic function 

I(k,co) may be expanded in a convergent Taylor series within a radius 6 about co = 0. 
(Note that this expansion is made possible by the analytic continuation of e + (k,co) in 
equation (21).) With co = pe*^, then 
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where 


I'CM) 


_ dI(k,co) 
doi 


Jw=0 


l"(k,0) = d I ( k > a) ) 


do;' 


a >=0 


p 277 o • fj 

The first integral is of the form \ e" which is identically zero. The next 

J ir 

term gives, in general, a zero imaginary part but a divergent real part; this term, how- 
ever, does not contribute to <|p(k)| 2 > because the imaginary part of I(k) is to be 
taken. In any case, the odd order derivatives of I(k,o>) would be zero for any even r 
(that is, any symmetric distribution function) since odd order derivatives will have odd 

integrands. The next term contributes 7 ri ^or, in general, ni the usual 

residue for a pole of order r - l'] and the remaining terms go to zero in the limit 6—0. 
For r = 4, the result is 




1 d Mk,a>)l = 47ri 


[e(M)] S 


dco 2 


f u>=0 


If the second moment of the distribution function is used to show that 
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^2 = 


rf^ 


r/i 


00 

where r(b) denotes the gamma function r(b) = \ dy y^ - -*- exp(-y), then 


P 2 (k) 


‘ > = ncr 2 
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/ON 

-N 
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\k 2 + A 
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> 

D 2 

r (i 


V 

V*/ 



(25) 


The concept of the Debye distance can be extended to nonequilibrium situations by 
finding a scale length in the equations whose use is analogous to D. In equation (25) the 
obvious choice is 


V T 

Do = 

2 20V, 


Yi\ 
\4 / 


M) 


= 1.48D 


20 



The quantity D 2 is later shown to be analogous to the Debye distance in that it serves 
as a scale factor for nonequilibrium screening of plasma fields. Equation (25) then 
simplifies to 


< 



> = 


no 2 



(26) 


As a final example, r = «> gives the square distribution fg(v) shown in figure 2, 
for which <|p(k,o))|^> can be integrated directly. If it is noted that 


f 


3 



for v < a and is zero elsewhere, then 



af 3 (v) 

8v 


^[fi(v + a) - 6(v - a)J 


If the dispersion relation 


is used, then 


so that 


o> k 2 = co p 2 + k 2 a 2 


e(k,co) 


= 1 + 


C0 T 


CO 2 - U) k 2 


Integration yields 

,2 


:|p 3 ( k )|' 


> = na 2 




1 + V. log (V^\ ♦ JV- 10 

kaa^k lu> k + kaJ 4kaco, 2 V^k " ka / 2a), 2 


(27) 


With O = — =- and an analogous "nonequilibrium Debye distance" D Q = — = \[3D, equa- 

0) k 6 U)p 

tion (27) simplifies to 
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2 

<|p 3 (k)| > = no 2 


(28) 



The Potential Energy Harmonics 

Figures 5 and 6 show experimental measurements of W(k), the Fourier harmonics 
of the total potential energy as a function of kD, a dimensionless wave number for two 
plasmas having initially Maxwellian |f(v) = f ^(v)j and square jf(v) = fg(v)j distribu- 
tions, respectively. Implicit in the comparison with theory is the fact that f j(v) and 
f 3 (v) are considered to be stationary distribution functions; this fact is restricted to one- 
dimensional plasmas as a result of the first-order kinetic theory considered later. 

Dawson (ref. 16), in fact, considers transient changes in nonequilibrium distribution func- 
tions, but application of these corrections is of limited usefulness in view of the time 
limitations of Dawson's treatment and the unnecessary complications incurred in using 
the modified distribution function. All plasmas described in this paper were generated 
with N = 1000 and AT = 0.1; the plasmas used to generate the data of figures 5 and 6 
both have nD = 20. This intermediate value of nD is large enough to allow nonequi- 
librium effects to be apparent and short enough to allow several modes to occur in the 
critical range of kD in a stationary state. Using kD, a dimensionless wave number, 
as the abscissa emphasizes the scale of separation between propagating, coherent oscilla- 
tions (kD < 1) and those disturbances damped out within a few Debye distances (kD > 1). 

In figure 5, the theoretical curve is obtained from equations (20b) and (24) and is 
shown as a continuous curve for clarity only. For a finite system size, the potential 

energy harmonics have meaning only for integral mode numbers c, where kD = - ■ 7rc ^ > . 

L 

The data for the Maxwellian plasma agree well with theory except for large fluctuations 
in the region kD < 1; these fluctuations are discussed below. The sum over the first 
20 modes shown here does, however, converge to within 10 percent of the predicted sum 
over these modes. 

For the plasma with f(v) = f^(v) (fig. 6), the theory curves for fg(v) and fg(v) 
calculated from equations (20b), (26), and (28) are both plotted. The term W 2 (k) is of 
interest not only because of its qualitative and quantitative similarity to Wg(k), combined 
with an analytically simpler form, but also because W 2 (k) may be more relevant on 
physical grounds. As will be shown in the section on relaxation effects, for moderate 
values of nD, transient relaxation of fg(v) results in a distribution close to fg(v) and 
the specific choice of W 2 (k) or W 3 (k) depends on scaling of first-order terms and 
averaging intervals. The scatter of the experimental data also makes any distinction 
meaningless. The agreement between experiment and theory is acceptable for short 
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wavelengths (kD > 1.6) but below this level not only are fluctuations large but also W(k) 
does not go to zero with kD; at best, there is some tendency for W(k) to peak at non- 
zero kD. 

There are several possible explanations of these difficulties, some based on lim- 
itations of the theory and others due to practicalities of the experiment. Several possi- 
bilities are as follows: 

(1) The energy balance of the long wavelength modes may be affected by nonlinear 
coupling not described by the linear theory; that is, the random phase approximation may 
not be applicable. 

(2) Long wavelengths may be excited by random (roundoff) errors in calculation of 
the particle positions. These errors are basic to the computer program used for these 
experiments which only approximated by the Ej term in equation (5b) the time depen- 
dency of the electric field during a time step; these errors could be eliminated only by 
using an exact program which is seriously limited by time considerations for nD > 10. 
Random position errors affect the long wavelengths preferentially, because with n(x) 
constant (that is, random) |p(k)| 2 = Constant; thus, 

2 

2 _ < | p (k) | > Constant 

* “ k 2 = k 2 

which diverges at small k. 

However they may be excited, these waves cannot be Landau damped, because, as 
pointed out by Weitzner (ref. 23), there are no particles with velocities equal to the phase 
velocity of (that is, synchronous with) the waves. Because of the limited number of parti- 
cles that can reasonably be treated in the computer experiment, this' cutoff is at least as 
high as kD * 0.3 and often higher. Also, because there is generally a small number of 
particles with high velocity, the poor statistics resulting from this scarcity upset the 
detailed balance of excitation and damping. 

(3) The zeroeth-order kinetic theory with which the calculations are compared is 
not applicable at long wavelengths for finite nD (nonzero graininess); that is, the limits 
g — 0, k — 0 do not commute. This problem has been investigated by Kaufman (ref. 24) 
who predicts a finite value of E^ 2 as kD — 0. The potential energy in a given model 
depends on a balance between excitation and damping to successive orders in the graini- 
ness parameter; that is, 

e„(k) + ge.(k) + g 2 e ? (k) + . . . 

W(k) oc 0 - - 1 2 

d 0 (k) + gd x (k) + g 2 d 2 (k) + . . . 

where e and d^ represent the excitation and damping coefficients to order y in g. 
In thermal equilibrium, this condition presents no difficulty because this ratio to each 
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order is 1/2, and only the fluctuations disturb the long wavelengths. For nonequilibrium 
plasmas, only the zeroeth- order contribution vanishes and other contributions will go to 
a finite level. 

The first-order contribution to the potential energy comes from collisional effects. 
By use of a phenomenological model for the collisions, deviation from the zeroeth- order 
theory is predicted for kD < (log e where X ~ nD. For nD = 10 and 60, the 

limits of the present experiments, the breakdown of the Vlasov theory comes at *0.66 
and 0.49, respectively; thereby the difficulty of adequately representing the long wave- 
lengths with a plasma of reasonable size (due to the L/D » 1 restriction) by the Vlasov 
theory is emphasized. Tripling the upper range of nD (to 180) and taking N = 3000 
(to maintain the same L/D) should make the Vlasov treatment good down to kD = 0.44. 
The improvement in resolution would hardly be justified in view of the greatly increased 
calculation time. Difficulties with the zeroeth-order (that is, Vlasov) results in no way 
prejudice the credibility of first-order results, for it is precisely the effects of nonzero 
graininess (that is, finite nD) that destroy the correlationless nature of the Vlasov equa- 
tion. It must be pointed out that this discussion puts the wave amplitude in doubt, but 
offers no insight into the large fluctuations observed. 

(4) Whatever the means of excitation of the long wavelength modes, Dawson (ref. 15) 
shows that they require progressively longer times to reach a stationary state as k — 0. 
Physically, Landau damping has been visualized as spontaneous emission of a plasma 
wave (plasmon) by a fast particle; the resultant energy loss is reflected in a slowing down 
of the particle (ref. 16). Overall balance is maintained by absorption of the wave by a 
particle whose velocity is nearly synchronous with the phase velocity of the wave. The 
scarcity of particles in the high-velocity tail of the distribution function means that sta- 
tistics on the long wavelengths will be poor; thus, the particular modes of the spectrum 
appear to be overexcited for long times. As expected, no major change can be seen in 
the total energy of the system during these fluctuations; energy is conserved, but is not 
distributed correctly. 


Electric-Field Spatial Correlations 

Theory.- By integrating over both co and k, it is possible to smooth the effect of 
long wavelength fluctuations and yet still retain enough information to distinguish between 
equilibrium and nonequilibrium results. The spatial correlation of the electric field, 
that is, 

C(X,0) = <E(x,t) E(x+X,t) > x ^ t 

results from such an integration. All details of interest are found in the first few Debye 
lengths, within the limits of applicability of the zeroeth-order theory as outlined by 
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Kaufman (ref. 24). This function is closely related to the Debye shielding in that the field 
at a point x + A can be shown to be correlated to the field at a point x, through an 
exponential- shielding law in which the Debye distance or its analog is the scale length. 

By using equation (20c), C(A,0) is easily evaluated for f^(v) and f2(v). For f j(v). 


Ci(A,0) 


L C °° D 2 exp(ikA) ~ 
47 tD 2 1 + k 2 D 2 


Closing the k-axis contour at infinity in the upper half plane gives, from the residue of 
the pole at k = +ikj^, 


c l( x,0) = i exp(- |) 


(29) 


Evaluation of C2PO involves a second-order pole. Once again, D2 = 1.48D is 
taken as a characteristic scaling distance, so that 

k ^ 00 k^Dg^ exp(ikA) 


C 2 (M) = 


4ttD 2 


(i + (kD 2 ) ! 

Evaluating the contour in the upper half plane yields 


dk 


(3 ° 

The quantity C2(A,0), which is quantitatively close to €3(^,0) is used for comparison 
of experimental results. As expected, integrating over k gives a smoother result than 
W(k), the potential harmonics, because correlations are averaged over more variables. 

Measurements .- Figures 7 and 8 show measurements of the variation of C(A,0) 
with A/D and A/D2 for f^(v) and fg(v), respectively, where nD = 20 for both 
plasmas. It should be noticed that for A = 0, C(A,0) reduces to potential energy, and 

thereby suggests the intepretation of Debye screening as the shielding of potential fields 
by the dielectric plasma medium. 

The behavior of C(A) for small values of T is erratic because of the highly 
artificial state in which the calculations are started. The plasma particles are initially 
distributed uniformly; electrons and ions alternate at equal separations, and thus only 
the shortest wavelength is initially excited. As the plasma is released and the proper 
correlations are established, longer wavelengths are excited (by graininess effects, non- 
linear interactions, and so forth) so that the W(k) spectrum fills in from the highest 
values of k, and gives for all distributions a transient appearance similar to Wg(k). 


25 



ilium i 


n in m mu mi ■■■miii 


I I I II IIUIIIM I ■■■ ■■■■■■ 


The resultant form of C(A,0) is shown in figure 7(a) for a Maxwellian plasma, together 
with the theoretical value C 1 (X,0) (given in eq. (29)). 

As correlations are established in the plasma, the W(k) spectrum and, therefore, 
C(A,0) reaches a stationary state characteristic of the average distribution function of 
the plasma and this state persists as long as the distribution function remains unchanged. 
For a Maxwellian plasma this stationary state corresponds to the expected Debye- Huckel 
shielding. Figure 7(b) shows the form of the shielding when the steady state has been 
reached and gives good agreement with theory; figure 7(c) shows a long-time average. 
Both spot checks and long-time averages give the correct result, and all times shown are 
less than the thermalization time for a nonequilibrium distribution. 

The experimental data for f(v) = f Q (v) can be compared with C 9 (X,0) both quali- 

0 i 1 2 * \ 9 

tatively and quantitatively not only because W 2 (k) closely approximates Wg(k) 

but also on the basis of relaxation effects treated in the next sections. A sequence of 
measurements similar to figure 7 is shown for the electric-field spatial correlation for 
f(v) = f 3 (v) in figure 8. The expected antishielding (that is, negative correlation) appears 
after moderate time (fig. 8(b)) and persists on a long-time average. Again, all times 
shown are less than the thermalization time. The agreement with C2(X,0) is, however, 
only qualitative. A major factor in the quantitative discrepancy is, once again, the over- 
excitation of long wavelengths evident in figure 5. Because of this overexcitation, the 
total potential energy (and, therefore, C 3 (A, 0 )) is too high, and thereby shifts the entire 
level of the curve. 

FIRST-ORDER THEORY AND MEASUREMENTS 
The Test Particle Theory 

The concept of the test particle was introduced by Rostoker and Rosenbluth (ref. 25) 
as a model for treating first-order (graininess) effects in a plasma. If one artificially 
tags individual particles in the plasma, the kinetics of these particles interacting with the 
field, or background, summed over the plasma describe the first-order behavior of the 
plasma. Such tagging of particles is impossible in a laboratory plasma, but the ease with 
which this can be done in a computer experiment illustrates the great versatility of this 
technique. 

The test particle is described by a distribution function f^(v) which is a subset of 

the overall one-particle distribution function; that is, f^(v) = f(v). In three dimensions 

i 

the Balescu-Lenard equation (ref. 17) is the kinetic equation for the one-particle distribu- 
tion function to first order in the graininess parameter. Eldridge and Feix (ref. 14) have 
solved the one-dimensional form of this equation starting from the first-order kinetic 
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equations (similar to eqs. (12) and (13)) written for the test particle distribution; that is, 


9fi(l) 


8t 


- T 1 < 7 i n J ""‘a *2 s s n f 1 - x 2) PijM < 31 > 


and 
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cr.u. 

9v« 


fid) fj ( 2 ) 


(32) 


where n^ is the normalized density n^yhg, where f(l) denotes the one-particle dis- 
tribution function of particle 1, and so forth. 

Collective effects are contained in the pair correlation function P(l,2), the form 
of which depends on the entire distribution function; because of this dependence, the test 
particle results emerge as the interaction of discrete particles with a continuum plasma. 
The two-particle interaction described by this pair function therefore emerges as the 
lowest order correction in a plasma of finite graininess. In equation (32), fj(2) is taken 
to represent the entire plasma with which the test particle (of type i) is interacting. 

By using the methods of Lenard (ref. 26), the Fourier transform of equation (32) 
can be substituted into equation (31); integration over k, the wave number, yields 



Metaequilibrium 

The most striking result of this one -dimensional kinetic equation follows from sum- 
ming over all test particle species. If both sides of equation (33) are multiplied by rq 
and summed over i, the two terms in braces cancel identically, from which it 
follows that 
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( 34 ) 


and there is no thermalization of the one-particle distribution function to first-order time 

r)T\ 

scale r 1 = — . This result holds only for a one-dimensional system, but is valid for both 

X Wp 

equilibrium and nonequilibrium plasmas and predicts a stationary state or "metaequilib- 


rium" for any initially stable distribution on a time scale r ^ ~ ^ « ^51. 


This pre- 


diction has served as the basis for checking the experiments with fg(v) against a non- 
equilibrium theory. If f j(l) is not summed over the entire distribution function, that is, 
f j(l) * f(v), changes in the test particle distribution do appear on a first-order time scale, 
and relaxation effects on the test particle distribution can be calculated and measured. 


These results bear some analogy to a one- dimensional gas of particles having only 
short-range interactions (binary collisions). In one dimension, conservation of energy 
and momentum dictates that the two partners in a binary collision retain their original 
velocities or at most swap velocities. If the test particle distribution is the entire one- 
particle distribution, no net change occurs; however, if only one of the colliding particles 
is a test particle and its collision partner is a "field" or background particle, changes in 
the test particle distribution will occur. 


In a plasma the interaction is collective (N-body) and makes this result somewhat 
surprising. First-order effects are, however, similar to the binary collision case in 
that energy is exchanged between particles by emission of a plasmon wake by test parti- 
cles, propagation of this wake through the plasma, and reabsorption by another particle 
of nearly the same velocity (ref. 16) . This emission and absorption gives a somewhat 
binary character to the interaction, whereas the propagation and damping (that is, absorp- 
tion) mechanisms are collective effects. This process is effectively a "collision in veloc- 
ity space, " and thus carries over the analogy of the hard sphere gas to the metaequilib- 
rium result described. Thermalization is expected to occur in a higher order descrip- 
tion of the plasma including ternary correlations - again in similarity to the short-range 
collision model where three -particle interactions would bring changes in f(v). The same 
result also follows directly from Balescu-Lenard equation. In three dimensions, changes 
in the distribution function do occur on a scale Tj by interaction of plasmons of wave 

vector k with test particles for which k • (vj - V2) = 0. The Fourier transform of 

the interaction potential is the same for one and three dimensions; thus, the equation may 

be one-dimensionalized by dropping the vector sign and changing a multiplicative constant 

(ref. 14). In one dimension, the wave vector is always parallel to particle velocities- 
3fj(v) 

thus, — = 0 if fi(v)=f(v). 
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The persistence of the metaequilibrium is shown in figure 9 for four plasmas with 
initial velocity distributions of the form f 3 (v) and values of nD of 10, 20, 40, and 60; 
all four velocity distributions are shown at a first-order time = 4. The dashed 

curve represents an equilibrium (that is, Maxwellian) distribution having the same energy 
and the corners marked at V = ±1.7 show the corners of the original square. In all 
cases, = 4 is much greater than the zeroeth-order time scale Wp'^; thus, qualita- 


tive differences between the plasmas at 
higher order effects. 


Tj are the result of secondary and possibly 




co T 


and the plasma has already substan- 


For nD = 10, r = 4 = 0 . 4 t 2 where T 2 = 

tially thermalized. In contrast to this, for nD = 60, r = 4 = O.O 6 T 2 (that is, <<T 2 ) 
and the original square shape is well retained. The changes that do occur, however, 
appear to be most rapid in the high velocity tail of the distribution, which tends to over- 
excite the long wavelength modes of potential energy; this overexcitation, to some extent, 
may explain the difficulty in obtaining quantitative agreement for |p(k)|^ and C(X). All 
distributions are modified from their original shape and, for example, for nD = 20, f(v) 
is very close to f£(v) which, as pointed out earlier, appears to be an intermediate form 
of f 3 (v) in its eventual thermalization. 


Relaxation of Test Particles 

If f i(l) is taken as a subgroup of the total distribution, for example, all particles 

af i (i) 

with velocities Vi such that V 1 -AVSV i ^V 1 + AV, * t ■ * 0 and equation (33) can 
be viewed as a Fokker- Planck equation of the form 

^ + -| F [x(v)Ii(l)-|-| r <'(v)fi(l)]=0 (35) 

A <Vi > A (<Vi 2 > - <Vi> 2 ) 

in which y(v) = — — — - is the drag and i^(v) = — ^ is the velocity diffu- 
At At 

sion coefficient. Such a test particle distribution is shown at initial time in figure 10(a) 
with its parent distribution f(v) (that is, the total velocity distribution at the same time). 
In this example, Vj = 0.4, AV = 0.05, and the parent distribution was initially a square 
[f(v) = f 3 (v); nD = 20j whose corners are shown at V = ±1.7. As time progresses, the 
test particle distribution interacts with f(v) and changes as shown in figures 10 (b) and 
10(c). The amplitudes of the time- averaged test particle distributions are amplified by 
factors of 15 and 25, respectively, in figures 10(b) and 10(c) to aid visual comparison with 
the parent distribution. 

Most of the particles seen by the test group have velocities < Vj; thus, the drag 
experienced by fi(v) makes <Vi> tend toward zero. This is true, however, only on 
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the average, and individual members of the group gain or lose more energy than the 
group average; this condition leads to dispersion or "diffusion" in velocity space. Both 
effects are indicated in figure 10(b). This dispersion in no way represents the effects of 
propagating errors in the calculations but is the expected microscopic behavior. After 
many plasma periods the test particles have relaxed to a stationary state which is not 
Maxwellian but a duplicate of (that is, synchronized with) the parent distribution. The 
slight discrepancy between f(v) and f^(l) in figure 10(c) is due to a slight shift in fj(v) 
arising from nonzero < Vi >; in this case <Vi> = 0.1. This result again supports the 
metaequilibrium concept in that the test group has not really "thermalized" (in the sense 
of becoming Maxwellian) but relaxes to a stationary nonequilibrium form. 

The exact form of the distribution function is more closely examined in figure 11. 

In figure 11(a) the parent velocity distribution is plotted with the theoretical curves for 
the square, Druyvesteyn, and Maxwellian distributions. (Once again, only the outline 
corners of the square are shown.) The departure from the initial square is evident, but 
the plasma has not yet become Maxwellian. The relatively good fit to the Druyvesteyn 
gives more quantitative justification for using the theory for f 2 (v) with the zeroeth- 
order experiments. Figure 11(b) compares the test particle distribution with the theory 
curve for f 2 (v). The test particle distribution fj(v) shows more fluctuations than 
f(v) but, once again, f 2 (v) provides a reasonable fit. 


Diffusion Coefficients 

Theory .- Quantitative checks of the first-order kinetic theory are, to some extent, 
better than those of the zeroeth-order theory. The diffusion of a test particle can be 
calculated in a very general way by following the work of Thompson and Hubbard 
(refs. 27 and 28). The equation of motion for the test particle is 

W = " H E iW 

where E^(t) is viewed as the stochastic field caused by the rest of the plasma as seen 
by particle i. For short times, the differential motion is given by 


Av i = " H J Q 1 E i\r (t " s )> 1 ‘ ®]' ds 


(36) 


and results in a diffusion coefficient 


^(v) = ^ Ai (v) = ^<Avi Avi> = 2 <vq Av A > = 2^J < 1 ds|E i (x,t) Ej[x(t - s), t - sjj: 


(37) 
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The time over which the correlations persist is short compared with the time during 
which a particle trajectory is significantly modified; thus, x(t - s) ~ x(t) - vs where v 
is the average particle velocity. Equation (37) then assumes the form of the autocorrela- 
tion coefficient C(A,r) with X = vs and t - s. The short correlation time also per- 
mits the upper limit of the s-integration to be extended to infinity, and, as a result, 


> exp i(kv + u>)s 




ds H(s) ~ < 


_ 1 _ 

277 


< 

|E(k,o>)| ' 


> exp i(kv + u))s 


where H(s) is the Heaviside unit step function. Performing the s-integration yields 

‘r. do) y dk <|E(k,o>)|^> 6(kv + u>) (38) 


In terms of the dimensionless velocity V = — and the first-order scaling time = — , 

I P 

one obtains a normalized diffusion coefficient 


^(V) = JL[<V 2 > - <V>5 = f do) f dk ^ 
^ dnL J v T J_ oo J -00 .3 i /,„ 


k J e + (k,cti) 


5(0) + kv) 


(39) 


The simplest case to treat is Vj = 0, for which the drag, a velocity-dependent force, 
goes to zero. (Because of the experimental necessity of nonzero AV around Vj, the 
dependence of drag on V introduces additional scattering not strictly due to diffusion. 
This effect could be more pronounced for Vj =£ 0.) The integration over io for Vi = 0 

then reduces to evaluation of <|E(k,w)j^> at a> = 0; that is, 


Mv) 


_ ^p 2 r a 
v T J-« 


_f(0) 


! ||e(k,0)j : 


dk = 


_ w p 


•I' 


: (k, 0) | ‘ 


dk 


(40) 


where f(0) = A. For f-^ = Aj exp^-j3jV^, 


J -°o o - v 


( 41 ) 


so that with A, = 


^ 277 Vrp 
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( 42 ) 


2a)p2 n co 


r"~ v^tv 1 2gj 2 

\ ^ dk = -JL = 0.4 

>Jn 


xb. (V) = £ .... 

\Z2irnv T 3 ^0 (l + k 2 D 2 ) 2 |k|k 2 fen 

Similarly, for f2(v) = A2 exp^-^v^ where A2 is defined as 


2y3 2 l/4 2 


r(i 


a 2 = 


and ^2 is defined preceding equation (25), then 

e (k,u>=0) = 1 + 


[ r (!3 


1/2 


'(!) 


aV2 


k^D2 2 


where D 2 = 1.48D. From this relation, 




ni/2 


^ 2 (V) - 2 



rf'l'] 
.4/ 


3/2 


= 0.687 


Finally, for F3 = — for ] v | < a = 3vrp, 

e(k,o>=0) = 


1 + ( k ^3, 


9 / a \2 

and D3 = — ) ; thus, 
6 \ w p/ 


^ 3 ( y ) = y = 0.866 


(43) 


(44) 


(45) 


(46) 


Measurem ents.- The experimental results of test particle diffusion calculations are 
shown in figures 12 and 13. By selecting the test particles in each plasma approximately 
three plasma periods after the initiation of the experiment (Tq = 20u>p-4j, the variance of 
the velocity is plotted as a function of time in units of the zeroeth- order time scale (that 
is Tq = Wpt^ for four plasmas having nD = 10, 20, 40, and 60 for an equilibrium plasma 
(fig. 12) and a nonequilibrium plasma, f(v) = f 2 (fig. 13). The initial slope of these 
lines gives the desired diffusion coefficient; in taking the slope, the relaxation time 
t W n t 

— = — t- must be used, rather than the zeroeth-order time scale on which the data are 
T 1 nD 
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plotted for clarity. The initial slopes given in table I for the equilibrium plasma are 
relatively independent of nD and check well with the theoretical value of 0.4. The non- 
equilibrium case, which in figure 11 showed increasing stability of the metaequilibrium 
with increasing nD, again shows the effect of nD on the thermalization time. For 
nD = 10, the substantial thermalization suggested in figure 11 is borne out by a diffusion 
coefficient of 0.5 which is close to the equilibrium value of 0.4. The relative stability 
of the initial square for nD = 60 (fig. 11) is reflected in a diffusion coefficient of 0.9 
compared with the theoretical value of 0.86 for fg(v). For nD = 20, the diffusion coeffi- 
cient of 0.65 indicates partial thermalization; the proximity of this value to the theoretical 
value of 0.70 for fg(v) again justifies the use of |E£(k)p and C^(A,0) theory for 
experiments with f 3 (v). The diffusion measurements are summarized in table I. 


Drag Coefficient Measurements 


Figure 14 shows drag measurements on the same four equilibrium plasmas studied. 
The average velocity of the three test groups /<Vi>q = 0.4, 1.0, and 1.4) is shown as a 

CUv^t 

function of normalized time the drag coefficient is given as the slope of the 

curves, and the straight lines represent the theoretical values calculated in reference 12 
from equation (33). The data generally agree with the theory except for the fluctuations 
in the nD = 40, = 1.0 case. Figure 15 shows on a much longer time scale the relax- 

ation of test particles initially at the thermal velocity. The extensive overlap of the data, 
even down to very small velocities, verifies the correctness of the relaxation time scale 
Tj, because the zeroeth-order times ^that is, used for the generation of the plasma 

vary by a factor of 6 as nD goes from 10 to 60. 


The data for f 3 (v) are more scattered, but generally exhibit the trend shown in fig- 
ure 16 for Vj = 1.0, where the straight line again represents the equilibrium drag. The 
thermalization of the plasma for which nD = 10 is again shown by its agreement with the 
equilibrium line. For nD = 20, 40, and 60, the drag (that is, slope of the data) is greater 
than that for the Maxwellian, as would be expected from equation (35). Figure 17 again 
shows the consistency of the long-time scaled behavior of the drag for the four plasmas 
considered. 


CONCLUSIONS 

In conclusion, an approximate, self-consistent N-body calculation has been used to 
study the microscopic behavior of one-dimensional plasmas in a controllable, physically 
realistic manner. The physical consistency (for example, energy conservation) in the 
approximate model improves as nD (particle density times Debye distance) increases 
because of smoothing of the graininess effects even though more crossings occur per 
fixed time step. A large range of values of the graininess parameter g has permitted 
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separation of collective and individual particle effects through the zeroeth- and first- 
order theory. 

By using the Vlasov theory, calculations of W(k), the Fourier harmonics of the poten- 
tial energy, give good agreement for equilibrium plasmas, but only fair agreement for 
nonequilibrium. The balance of long wavelength modes is probably not well described by 
the Vlasov theory for finite g and in this range, even the equilibrium agreement may be 
fortuitous. Integrating over wavelength, the electric field spatial correlation C(A) shows 
excellent agreement in equilibrium and qualitative agreement out of equilibrium; in par- 
ticular, antishielding is observed as a stationary effect for a range of several Debye 
lengths. 

The metaequilibrium or stationarity of nonequilibrium distributions was verified and 
serves as the basis for study of nonequilibrium plasmas. The test particle approach is 
used for calculation of first-order effects and produces good agreement with the expected 
drag and velocity diffusion. The diffusion coefficient is shown to be a good indication of 
the extent of thermalization in a nonequilibrium plasma. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Langley Station, Hampton, Va., June 13, 1968, 

129-02-01-01-23. 


34 



APPENDIX 


GENERAL DESCRIPTION OF THE COMPUTER PROGRAMS 

The data in this report have been generated by three separate computer programs. 
Program I generates the time history of the particle motions in a one-dimensional, self- 
consistent electrostatic plasma. Program II uses the position data from program I to 
calculate the Fourier modes of the potential energy and the Debye screening. Pro- 
gram III uses the velocity data from program I to give the total and test particle distribu- 
tion functions and the drag and velocity diffusion coefficients. 

The data from program I, in the form of position, velocity, and particle order arrays 
for successive time steps (that is, cycles) are output on magnetic tape to facilitate their 
later use with programs n and III. Kinetic and potential energy and total momentum are 
also output at each cycle and serve as checks on self-consistency and disaster monitors. 
By using three separate programs, the diagnostic parameters of programs II and III may 
readily be changed without repeating the relatively time-consuming trajectory calculation. 

In program I, equations (5a) and (5b) are used to calculate the positions and veloc- 
ities of identified particles in a stepwise manner; that is, the xj(Tj + AT) of one cycle 
becomes the x^(T) for the next cycle. The electric field is calculated using equa- 
tion (2b) (a sum of N + 1 terms) rather than equation (2a) (a sum of N(N + 1) terms) 
to compute the field seen by each particle. The sole disadvantage of this method is the 
required generation and storage of a particle-order array in which the identification num- 
ber of all the particles are stored in order of increasing position value from x = 0 to 
x = L. At the end of each cycle, the particle order is established by using a subroutine 
in which the order from the previous cycle is used as an approximation to the new order. 
If AT « 1, few crossings will have taken place and only minor rearrangement is 
necessary. 

The initial particle velocities should be assigned with a random number generator to 
avoid any correlation between position and velocity, which would be unrealistic. The ini- 
tial positions should, however, be assigned in a regular manner to avoid large excesses 
of potential energy. As pointed out in the discussion on potential energy harmonics, a 
completely random distribution of charge gives far too much potential energy at long 
wavelengths. This situation increases the time necessary for proper correlations to be 
established in the plasma because long wavelengths are weakly damped. Regular initial 
ordering of the particles provides precise control on the initial potential energy which 
may be concentrated on any wavelength or combination thereof, rather than on the short- 
est wavelength, as was done in these experiments. The transfer of potential energy from 
one wavelength to another is an interesting problem in itself. 

Program II uses the tape output from program I as input to calculate the Fourier 
modes of potential energy and the spatial correlation of the electric field. The particles 


35 



APPENDIX 


(that is, sheets) are taken to be infinitely thin so that the symmetric Fourier transform 
of the charge density at time T is given as 

p<k ’ T) ■ f C <fe elkX 2 ff i 6 [ x - x i (T 3 - H 1 °i eikXj<T) 

j=l j =1 

where Oj and Xj(T) denote the charge and position (at time T) of the jth particle and 

2 7TC 

the wave number is k = — — . The total potential energy modes are, therefore, 


t 2 N 

w(k) = i52 I I 0j[cos2kXj(T) + sin 2 kXj(T)| 

* 1 T=Tj_ j=l 

where AT' is the time interval between samples of W(k). The extent of the time aver- 
age (that is, T 2 - Tj) depends on the wave numbers being investigated because the 
Landau damping goes to zero with k; thus long wavelengths take longer to come to equi- 

■pv 

librium. A minimum averaging time of Ti = — — was used in all calculations and one 

1 w p 

needs no more than one calculation of W(k) in a time l/u>p if redundancy is to be 
avoided. Accuracy at long wavelengths is still limited, however, by required averaging 
times which can become of the order of T 2 , in which time the nonequilibrium spectrum 
changes, and thereby voids the calculation. Also, random position errors inherent in an 
approximate program continually lead to anomalous excitation of long wavelengths, and 
make the problem of obtaining good statistics on long wavelengths even more difficult. 


The calculation of the electric-field spatial correlation is straightforward. By using 
the particle order data output by program I, the electric field is calculated at regular 
intervals from x = 0 to x = L. If the basic interval is l and X is some multiple 
of l, then 

^2 n max 

C(A,0) = ~^~2 X X E(nZ,T) E(nZ + A,T) 

2 1 2D T=Ti n=Q 


where n max is the integer part of (L - X/l). A slightly modified version of this pro- 
gram has been used by Massel (ref. 29) to calculate the spatial and temporal correlations 
of potential fluctuations. 

Program III provides the most sensitive diagnostics on the approach to equilibrium 
in the one- dimensional plasma. The overall distribution function is calculated by sorting 
the velocity arrays. The problem of time averaging is, however, a delicate one and care 
must be taken in interpreting the error (that is, standard deviation) in the measurements. 
The first-order relaxation time Tj defines the time for statistical independence of mea- 
surements of f(v). If, for example, n independent measurements each yield mj 
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counts in a given velocity interval, then the probable error in the measurements is 
m = \i y m^. Thus measurements of f(v) should be separated by times >tj. 


Because the distribution function thermalizes in a time one must also require 
that all components of an average be taken in a time « Tg. It is often difficult to satisfy 
both requirements; thus several measurements must be taken within a time r^, at the 
expense of incurring a larger measurement error because the samples are not indepen- 
dent. At the very worst, if all the samples were completely dependent, the error would 


be VnN(v) = 


n 


' n 2 y 

i=l 


irq = nl 



The actual error is between these two values, but 


for small nD plasmas, one can expect rather marginal improvement from time 
averaging. 


Measurement of the diffusion and drag coefficients is made with the same sorting and 
averaging technique as used for the total distribution function. By selecting particles 
within a narrow velocity interval, the mean and variance of the test particle distribution 
as a function of time yields the drag and diffusion. If groups in the same velocity range 
are selected at different times T, the averaged test particle distribution as a function 
of r (the time elapsed since selection) is given as 


fi(V,r) = 



fi(V,T + r) 


For a test particle calculation, however, statistical independence is achieved more 
rapidly because the particles are identified. To see why independence is achieved more 
rapidly, consider the changes produced by a one-dimensional binary collision in which the 
two particles either swap velocities or reverse the signs of their velocities, because of 
conservation requirements. Such an interaction changes nothing in the overall distribu- 
tion function because only the quantity of particles in a velocity interval is noted. If the 
particles are identified, as in a test-particle calculation, the state of the sample has 
changed and some degree of independence is attained. Because the number of test parti- 
cles is a small fraction of the field particles, with which almost all interactions occur, the 
time for statistical independence of samples should probably be reduced by a factor 
depending on the ratio of field particles to test particles. Thus, the expected accuracy of 
a test-particle calculation can be very high, even with closely spaced samples. 


37 


REFERENCES 


1. Alder, B. J.; and Wainwright, T. E.: Studies in Molecular Dynamics. 

I. General Method. J. Chem. Phys., vol. 31, no. 2, Aug. 1959, pp. 459-466. 

II. Behavior of a Small Number of Elastic Spheres. J. Chem. Phys., vol. 33, no. 5, 

Nov. 1960, pp. 1439-1451. 

2. Lomax, R. J.: Transient Space- Charge Flow. J. Electronics Control, vol. 9, no. 2, 

Aug. 1960, pp. 127-146. 

3. Birdsall, Charles K.; and Bridges, William B.: Space-Charge Instabilities in Electron 

Diodes and Plasma Converters. J. Appl. Phys., vol. 32, no. 12, Dec. 1961, 

pp. 2611-2618. 

4. Bridges, William B.; and Birdsall, Charles K.: Space- Charge Instabilities in Electron 

Diodes. II. J. Appl. Phys., vol. 34, no. 10, Oct. 1963, pp. 2946-2955. 

5. Burger, Peter: The Opposite- Stream Plasma Diode. SEL-64-012 (NASA Grant 

NsG 299-63), Stanford Univ., Apr. 1964. 

6. Buneman, O.: Dissipation of Currents in Ionized Media. Phys. Rev., Second ser., 

vol. 115, no. 3, Aug. 1, 1959, pp. 503-517. 

7. Goldstein, Charles M.: Electron Flow in Low-Density Argon Gas Including Space- 

Charge and Elastic Collisions. NASA TN D-4087, 1967. 

8. Hasegawa, Akira; and Birdsall, Charles K.: Sheet- Current Plasma Model for Ion- 

Cyclotron Waves. Phys. Fluids, vol. 7, no. 10, Oct. 1964, pp. 1590-1600. 

9. Dawson, J. M.: Investigations of Nonlinear Behavior in One- Dimensional Plasma 

Model. Symposium on Computer Simulation of Plasma and Many-Body Problems, 
NASA SP-153, 1967, pp. 25-30. 

10. Buneman, O.; and Kooyers, G.: Computer Simulation of the Electron Mixing Mechan- 

ism in Ion Propulsion. AIAA J., vol. 1, no. 11, Nov. 1963, pp. 2525-2528. 

11. Yu, S. P.; Kooyers, G. P.; and Buneman, O.: Time- Dependent Computer Analysis of 

Electron- Wave Interaction in Crossed Fields. J. Appl. Phys., vol. 36, no. 8, Aug. 
1965, pp. 2550-2559. 

12. Anon.: Symposium on Computer Simulation of Plasma and Many-Body Problems. 

NASA SP-153, 1967. 

13. Eldridge, O. C.; and Feix, M.: One-Dimensional Plasma Model at Thermodynamic 

Equilibrium. Phys. Fluids, vol. 5, no. 9, Sept. 1962, pp. 1076-1080. 

14. Eldridge, O. C.; and Feix, M.: Numerical Experiments With a Plasma Model. Phys. 

Fluids, vol. 6, no. 3, Mar. 1963, pp, 398-406. 

38 


I 



15. Dawson, John: One- Dimensional Plasma Model. Phys. Fluids, vol. 5, no. 4, Apr. 

1962, pp. 445-459. 

16. Dawson, John M.: Thermal Relaxation in a One-Species, One- Dimensional Plasma. 

Phys. Fluids, vol. 7, no. 3, Mar. 1964, pp. 419-425. 

17. Montgomery, D. C.; and Tidman, D. A.: Plasma Kinetic Theory. McGraw-Hill Book 

Co., Inc., c.1964. 

18. Smith, Craig; and Dawson, John: Some Computer Experiments With a One- 

Dimensional Plasma Model. MATT-151 (Contract No. AT(3Q-1)-1238), U.S. At. 
Energy Com., Jan. 1963. 

19. Landau, L.: On the Vibrations of the Electronic Plasma. J. Phys. (U.S.S.R.), vol. 10, 

1946, p. 25. 

20. Jackson, J. D.: Longitudinal Plasma Oscillations. J. Nucl. Energy: Pt. C, vol. 1, 

no. 4, July 1960, pp. 171-189. 

21. Rostoker, Norman: Fluctuations of a Plasma (I). Nucl. Fusion, vol. 1, no. 2, Mar. 

1961, pp. 101-120. 

22. Feix, Marc R.; and von Hagenow, Karl: Connection Between Correlations and 

Fluctuations in a Plasma. Phys. Fluids (Res. Notes), vol. 8, no. 8, Aug. 1965, 
pp. 1565-1567. 

23. Weitzner, Harold: Plasma Oscillations and Landau Damping. Phys. Fluids, vol. 6, 

no. 8, Aug. 1963, pp. 1123-1127. 

24. Kaufman, Allan N.: Extension of the Plasma Kinetic Equation to Small Wave Numbers. 

Phys. Rev. Letters, vol. 17, no. 22, Nov. 28, 1966, pp. 1127-1129. 

25. Rostoker, Norman; and Rosenbluth, M. N.: Test Particles in a Completely Ionized 

Plasma. Phys. Fluids, vol. 3, no. 1, Jan. -Feb. 1960, pp. 1-14. 

26. Lenard, Andrew: On Bogoliubov's Kinetic Equation for a Spatially Homogeneous 

Plasma. Ann. Phys., vol. 10, no. 3, July 1960, pp. 390-400. 

27. Thompson, W. B.; and Hubbard, J.: Long-Range Forces and the Diffusion Coefficients 

of a Plasma. Rev. Modern Phys., vol. 32, no. 4, Oct. 1960, pp. 714-718. 

28. Hubbard, J.: The Friction and Diffusion Coefficients of the Fokker-Planck Equation 

in a Plasma. Proc. Roy. Soc. (London), ser. A, vol. 260, no. 1300, Feb. 7, 1961, 
pp. 114-126. 

29. Massel, Gary Alan: The Kinetic Theory of Forced Oscillations in Plasmas. Ph. D. 

Thesis, North Carolina State Univ., 1967. 


39 



TABLE I.- DIFFUSION COEFFICIENT FOR ZERO VELOCITY TEST PARTICLES 
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nD 

Theory 

fl(v) 

Experiment 

fj(v) 
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0.40 
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20 

.40 

.42 

.86 
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.65 

40 

.40 
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.86 
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60 

.40 
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.86 
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Figure 5.- Potential energy spectrum for f(v) = fi(v). N = 1000; nD = 20; Time averaged from 75 - BSttfn" 1 * 








(C) fj(V) at To < T < To + 


■!; <Vj> = 0.1. 


Figure 10.- Relaxation of test particles. 
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(a) Parent distribution f(V). 



(b) Test particle distribution fj(V) x 23. 

Figure 11.- Comparison of test particle and parent distributions with theory at T = 230cj p _1 . 
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Figure 12.- Test particle diffusion for f(v) = f^v). T 0 = ZOwp" 1 ; (Vj) 0 = 0. 




Figure 13.- Test particle diffusion for f(v) = T 0 = 
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Figure 14.- Test particle drag. f(v) = fi(v). T 0 = 20w D _1 . 
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Figure 17.- Test particle relaxation. f(v) 
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